Mean-field theory of a plastic network of integrate-and-fire neurons 
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We consider a noise driven network of integrate-and-fire neurons. The network evolves as re- 
sult of the activities of the neurons following spike-timing-dependent plasticity rules. We apply 
a self-consistent mean-field theory to the system to obtain the mean activity level for the system 
as a function of the mean synaptic weight, which predicts a first-order transition and hysteresis 
between a noise-dominated regime and a regime of persistent neural activity. Assuming Poisson 
firing statistics for the neurons, the plasticity dynamics of a synapse under the influence of the 
mean-field environment can be mapped to the dynamics of an asymmetric random walk in synaptic- 
weight space. Using a master-equation for small steps, we predict a narrow distribution of synaptic 
weights that scales with the square root of the plasticity rate for the stationary state of the system 
given plausible physiological parameter values describing neural transmission and plasticity. The 
dependence of the distribution on the synaptic weight of the mean-field environment allows us to 
determine the mean synaptic weight self-consistently. The effect of fluctuations in the total synaptic 
conductance and plasticity step sizes are also considered. Such fluctuations result in a smoothing 
of the first-order transition for low number of afferent synapses per neuron and a broadening of the 
synaptic weight distribution, respectively. 
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I. INTRODUCTION 

The brains of Hving animals are perhaps the most 
complex organs one can find. These are networks of 
neurons formed through the connecting synapses, and 
their proper functioning is crucial to survival. Theo- 
retical studies of the dynamics of neural networks have 
contributed to our understanding of how these networks 
might function P, 0]. Recent progress in the study of 
synaptic plasticity is opening up new opportunities for 
understanding how these networks can form [J, [3, 111- 
One of the most important aspects of a functional neu- 
ral network is that the strength of its connections can 
change in response to the history of its activities, that 
is, it can learn from experience 0, H B @]- Among the 
theories of neural learning, the best known, Hebbian the- 
ory, states that when a neuron partakes in the firing of 
another neuron, its ability of doing so increases An 
oversimplified version of the theory states that "cells that 
fire together, wire together" which is generally adequate 
for rate-based views of neurons [llj. However, the ac- 
tual activities of neurons are described by "spikes" , that 
is, action potentials that are produced when the state 
of a neuron meets certain criteria, for example, when its 
membrane potential reaches a threshold value 'l2!|. It was 
discovered that the exact timing of these spikes plays a 
great role in determining the synaptic plasticity, adding 
an element of causation to the correlation requirement 
of the learning process. That is, the presynaptic neu- 
ron must fire before the postsynaptic neuron for the for- 
mer to actually take part in firing the later. In fact, it 
was discovered that a sharp change from potentiation to 
depression of the synaptic efficacy can occur as the fir- 
ing time of the postsynaptic neuron precedes that of the 
presynaptic one [1, @, [U, [i3l ■ 

An important feature of Hebbian learning is the pres- 
ence of positive feedback. That is, the stronger a synapse 



is, the more likely that it will be potentiated. Often 
in simulated plastic networks, this presents a possibil- 
ity of runaway synaptic weights, which are often curbed 
through the introduction of cutoffs or other devices. 
In some schemes, the potentiation of a synapse is sus- 
pended when its weight exceeds a preset cutoff value 
til HE [II- While such cutoffs might be justified by 
physiological limitations of the cells, in simulation stud- 
ies they often result in a pile-up of synaptic weights dis- 
tributed near the cutoffs that is not observed in real bio- 
logical neural networks [H, [H, [H [13, [HJ, . Runaways 
or pile-ups can often be mitigated by softening the cut- 
off [m, [ij, [2^ . Alternatively, it has been proposed from 
experimental observations that the potentiation and de- 
pression processes could be asymmetric, such that po- 
tentiation is an additive process to the current synaptic 
weight while depression is a multiplicative process. Such 
a mechanism was shown to produce a uni-modal distribu- 
tion of synaptic weights free of pile-ups among the many 
afferent synapses of a single neuron when the inputs are 
driven by Poisson spiking sources [22l |. 

Besides viewing brains microscopically as networks of 
neurons as discussed above, efforts aimed at understand- 
ing the structure and function of a living brain also in- 
clude studies at the macroscopic scale, for example of 
the interactions between different anatomical regions or 
even of the role played by the entirety of the brain as 
a vital organ [10, [23|- The contrast of the two scales 
is very analogous to the study of condensed matter sys- 
tems, which has faced similar challenges of bridging our 
understanding of the microscopic with our observation 
of the macroscopic. Among various tools employed by 
theorists studying condensed matter systems, mean-field 
theory has been proven valuable in that it often allows 
one to obtain a quick grasp of what macroscopic states of 
the systems can be expected from the set of microscopic 
mechanisms they follow '^sj. While it is well understood 
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that mean-field theories could fail to predict correct scal- 
ing behavior of systems in critical states, where fluctu- 
ations and long-range correlations are important, they 
generally are adequate (and often the first step in a sys- 
tematic procedure) in qualitatively describing the system 
in stable phases, where fluctuations are of limited range, 
and useful in revealing the structure of the phase spaces 
of the systems. The later is especially desirable for biolog- 
ical systems where large numbers of empirical parameters 
and, consequently, vast phase spaces are often involved 
in microscopic models of the systems. 

In the current study, we consider a network of 
integrate-and-fire neurons [265 driven by Poisson noise of 
fixed frequency for all neurons. The interaction between 
the neurons is taken to follow the neural transmission 
model proposed by Tsodyks, Uziel, and Markram [23| 
(the TUM model), which can account for the saturation 
effect of neural transmitter and the short-term depression 
of synaptic conductance. While our approach can han- 
dle alternative choices, the synaptic weights between the 
neurons are allowed to evolve following the spike-timing- 
dependent plasticity rules proposed by Bi and coworkers 
[22 . [2^ . A mean- field theory is used to determine the 
self-consistent average firing rate of the neurons. The 
noise driven firings dominate the small synaptic weight 
regime, while self-sustaining firing activity is triggered in 
the large synaptic weight regime. As the mean synap- 
tic weight of the network is varied, the mean-field theory 
predicts a hysteresis for the transition between regimes of 
noise dominance and self-sustaining activity. Within the 
mean- field environment, assuming the neurons are firing 
with Poisson statistics, the dynamics of a single synapse 
can be viewed as a random walk process in synaptic- 
weight space. Under the small-jump assumption, we 
use the master equation to calculate the drift and dif- 
fusion coefficient for the random walk. The resulting 
Fokker-Planck equation allows us to predict the station- 
ary synaptic- weight distribution of the process and close 
the self-consistency with the requirement that the sta- 
tionary distribution reproduces the mean synaptic weight 
characterizing the mean-field environment. For the plas- 
ticity rules and ranges of the parameters we considered, 
the synaptic weights form a narrow distribution, hav- 
ing a width proportional to square root of the plasticity 
rate. Finally, we extend the mean-field approximation 
to include the effect of fluctuations in the synaptic con- 
ductance as well as the variations of jump sizes in the 
random walk of synaptic weights. 

Of course, despite our mention of brains of living crea- 
tures in our introduction, the current analysis cannot 
hope to address issues of dynamics involved in such a 
venue for myriad reasons. Among others, here we imag- 
ine homogeneous, stationary networks, which is surely 
not the case in brain. However, there are continuing, 
revealing experiments on cultured neural networks with 
perhaps several hundred individual neurons in which all- 
to-all coupling is not an unreasonable approximation or 
starting point [H, [H HO, [M HI . This class of exper- 



iments represents an important step and will continue 
to provide important insights into the behavior of more 
complicated networks. Our aims are to improve the un- 
derstanding of the stationary, statistical properties of 
such plastic network and ultimately address dynamical 
behavior during formation. Furthermore, cultured net- 
works are on a scale approachable by the modeling, nu- 
merical and analytic work such as presented here. By 
analogy with a variety of familiar statistical mechanical 
models, sufficiently large homogeneous networks in which 
the number of "neighbors" of a particular element grows 
proportionally to the number of elements in the ther- 
modynamic limit, are expected to be well-described by 
the type of mean-field analysis employed in this study 
(see, e.g., chapter 3 of [sS]). In a separate publication 
we will present results of extensive numerical simulations 
on integrate-and-fire (and other representative neuron 
model) networks, the results of which can be put in the 
proper perspective via comparison with the calculations 
presented here and with results of in vitro experiments. 

The remainder of this paper is laid out as follows. 
In Section [Til we describe the model of the plastic net- 
work, which consists of (i) the dynamics of the mem- 
brane potential, (ii) the model of synaptic transmission, 
and (in) the plasticity rules. In Section IIIIl we apply 
the mean-field method to a network with fixed synaptic 
weights and determine the state of the system through 
response functions deduced from the dynamics of neuron 
and synaptic transmission under a given mean synaptic 
weight. In Section llVl we map synaptic plasticity to a 
random walk process in synaptic-weight space, determine 
the stationary synaptic-weight distribution, and close the 
self-consistency in the mean synaptic weight of the net- 
work. In Section |Vl we expand the mean-field theory to 
include the consideration of fluctuations in total synap- 
tic current of a neuron and fluctuations in step sizes of 
synaptic weight changes. Finally, we summarize and con- 
clude in Section IVll 



II. MODEL 

We consider a noise-driven plastic network of integrate- 
and-fire neurons. The neurons are cou pled using the 
TUM model of neural transmission [27l. |34|. described 
below. The noise is modeled by randomly forced firing 
of the neurons following Poisson statistics at a fixed fre- 
quency. 



A. Integrate-and-fire neuron 

The integrate-and-fire model is a single-compartment 
neuron model where the state of a neuron i is described 
by a membrane potential Vi . The dynamics of the mem- 
brane potential follows the differential equation of a leaky 
integrator [351] 
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(1) 



where Tm is the leak time for the membrane charge, which 
is given by the product of the total membrane capaci- 
tance Cm and resistance Rm, while Vq is the resting po- 
tential when the neuron is in the quiescent state. The 
total synaptic current /sy„ is a sum over each afferent 
synapse j for the neuron i 



Isyn — 



R ^ 



Wj,iYj,i (Uj,i - Vi) , 



(2) 
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where, for the synapse connecting neuron j to i, 
the synaptic weight, TZj^i is the reversal potential for the 
ion channels, and Yj^i is the fraction of active transmit- 
ters. The dimensionless synaptic weight Wj^i can be inter- 
preted as the maximum synaptic conductance, achieved 
when Yj i — 1, measured in units of the membrane con- 
ductance of the neuron. In the current study, we 
consider cases in which there is only one type of ion 
channel at all synapses so that they share the same re- 
versal potential TZj.i = TZ. In these cases, the difference 
of membrane potential from the reversal potential can be 
factored out and Eq. ([1]) becomes 



dt 



where 



(3) 



(4) 



is the total synaptic conductance (in units of R^^). 

In this model, a neuron fires when its membrane poten- 
tial reaches a threshold value, Vth. Then, its membrane 
potential drops immediately to a reset value, VJ-. The ac- 
tion potential of the integrate-and-fire model is assumed 
to be instantaneous and is not modeled explicitly. The 
spike train produced by the neuron i is defined as the 
function 



Si = '^6 {t - ti^n) 



(5) 



where ti^n is the time when the neuron i fires for the n-th 
time. 



B. Tsodyks— Uziel— Markram model of neural 
transmission 

The fractions Yj^i of the active transmitters are de- 
scribed by the TUM model [l^] of neural transmission, 
where the transmitters are distributed in three states: 
"active", with the fraction Y; "inactive", with the frac- 
tion Z; and "ready-to-release", with the fraction X. For 



a synapse with presynaptic neuron j and postsynaptic 
neuron i, these fractions follow the dynamics [13] 
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where td is the decay time of active transmitters to the 
inactive state, tr is the recovery time for the inactive 
transmitters to the ready-to-release state, and u is the 
fraction of ready-to-release transmitters that is released 
to the active state by each presynaptic spike. With the 
conservation rule 



X, 



Y. 



(7) 



there are two independent variables per synapse. Since 
the multiplying factor of the spike train Sj in the dynam- 
ics of Xj^i depends on Xj^i itself, which is discontinuous 
when there is a (5-function in Sj due to a spike at given 
time i, we must specify how the value of Xj^i should be 
evaluated at the time of the spike. Consistent with the 
TUM dynamics, the values of the factors multiplying Sj 
at the discontinuities are to be evaluated immediately 
before the discontinuities. 

The TUM model is very flexible. Since the variable 
Y represents a fraction in the TUM model, its value can 
saturate. As a result, an increase in the firing of the 
presynaptic neuron can only produce a less-than-linear 
increase of the active transmitters. However, a linearity 
of active transmitters on the firing rate can be recovered 
in the u — > limit, maintaining the average level of the 
synaptic conductance by keeping the product Wj^u con- 
stant. Additionally, the presence of the inactive state, Z, 
mimics the short-term depression of neural transmission, 
where repeated firings of the presynaptic neuron over a 
short period of time will reduce the maximum value at- 
tainable by Y temporarily through depositing the trans- 
mitters into the inactive state before their recovery. Fi- 
nally, manipulating tr/td can reproduce a model with a 
single, saturating "species" Y (see below). 

Since the dynamics ^ depend only on the spike train 
of the presynaptic neuron j, without external disturbance 
all the efferent synapses of a neuron should have the same 
values of transmitter fractions. Thus, we can drop the 
subscript i from dynamics ^ and regard these fractions 
as the properties of the presynaptic neuron j . Such sim- 
plification is not applicable for transmitter dynamics that 
depend, for example, on the postsynaptic neuron, that 
have synaptic-weight-dependent parameters, or synapse- 
dependent noise. These complications are not considered 
here. 
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C. Noise 

While the integrate-and-fire and TUM dynamics are 
both deterministic, we model the stochasticity of the net- 
work with additional noise driven firing events following 
Poisson statistics with the frequency Xn for each neu- 
ron. The noise driven firings are treated the same way 
as threshold firings, that is, the membrane potentials are 
brought instantaneously to the reset value and the 
firing times are included in the spike trains ([5]) of the 
transmitter dynamics ([6|). 



D. Spike-timing-dependent plasticity 

It has been observed experimentally that the change of 
synaptic efficacy depends on the precise timing between 
the presynaptic and postsynaptic spikes [l3|. When a 
presynaptic spike precedes a post-syna ptic spike, follow- 
ing van Rossum, Bi, and Turrigiano [22, we take the 
synapse to be potentiated by an amount proportional to 
g-At/TA -ftrhere At is the timing difference between the 
two spikes and r4 is the size of the potentiation time 
window. Similarly when a postsynaptic spike precedes a 
presynaptic spike, the synapse is taken to be depressed by 
an amount proportional to e~'^*/'^B -^yjiere tb is the size 
of the depression time window [13, . One reasonable 
mechanism for the cell to determine the interval between 
spikes is to imagine the cell evaluates the concentration 
of a decaying chemical species released at the first spike 
[36| . Following such a mechanism, the synaptic changes 
proposed in 22] can be produced for the synaptic weight 



by the dynamical equation 



dw 



dt ^ 



i Bi Sj , 



(8) 



where A is the potentiation constant for causal spike 
pairs, r is the depression factor for anti-causal spike pairs, 
while ^s and Bs are the amplitudes of such (assumed) 
chemical species controlling the magnitude of potentia- 
tion and depression. These amplitudes are assumed to 
follow the dynamics [s^ 



da^ 
dt 



faS-j, a ^ A,B 



(9) 



where different choices of the spike increments /o- lead to 
different schemes of spike pairing. For example, /ct = 1 
leads to "all-to-all" pairing of presynaptic and postsynap- 
tic spikes while /o- = 1 — cr leads to "nearest-neighbor" 
pairing of spikes [1^ . We note the stochastic differential 
equations in Eqs. ([8]) and ([9]) above should follow Ito's 
interpretation |37i] . That is, when the value of a jumps 
because of a spike in S, the value of fa, which determines 
the size of the jump, is evaluated immediately before the 
jump. The effects of the two choices of are identical 
when spike pairs are far and apart but differ when pos- 
sible spike pairs are close compared with the decay time 



To- of the respective timing chemicals. With the choice 
fa — \, the values of ^ or B can increase linearly with the 
frequency of spikes without limit. On the other hand, for 
the choice = f — cr, the value of a will never exceed I, 
its value immediately after a single spike. It is reasonable 
to expect, in a realistic setting, that the values of A or i? 
will increase with increasing spike rates, but they should 
saturate and be limited to some physiologically deter- 
mined maximum values with further increase in spike 
rates. Assuming spike increments = '^o- (1 ~ c) with 
Ua < 1 can produce this behavior. The parameter Ua can 
be interpreted as the mean fraction of chemicals released 
at each presynaptic firing relative to the amount that can 
possibly be produced without a exceeding its maximum 
value of 1. 

Additionally, certain effects of fatigue can also be ex- 
pected in synaptic plasticity, [s^ That is, while a {—A 
or B) can saturate when the spike rate increases, the sat- 
uration value itself can be expected to decrease when a 
high spik ing rate is sustained for an extended period of 
time [271, |3J| . Such effect can be modeled by assuming 
an additional inactive state , which takes up a fraction 
of a that can be produced for each spike, so that the in- 
crement becomes fa = {1 — a — !„). One can assume 
simple dynamics for !„, given by 



dt 



G 
Tier 



A. 

TRa 



(10) 



so that Icr is fed by the presence of a with the rate rf^ 
and decays with a recovery time constant tro-- One ex- 
pects that the rate of a's feeding into the inactive state 
[first term on the right of Eq. PT!]) ] should be less than the 
rate of its own decay [first term on the right of Eq. ([9])], 
so we require the condition tjo- > t^. 

The considerations outlined above conspire to produce 
a variant of the TUM model similar to what we have used 
for the neural transmitters apart from possibly different 
values of the constants involved. In the current simula- 
tion study of plastic integrate-and-fire networks, we re- 
tain the qualitative structure of the system but simplify 
by assuming Ta = Tja = tu, tr^ = tr, and Ua = u so 
that the same values of Yi calculated for the neural trans- 
mitter in the TUM model can be used as surrogates for 
Ai and Bi in the synaptic plasticity yielding 



dw 



dt ^ 



ill) 



These simplifications reduce the amount of computa- 
tion required without sacrificing qualitative and semi- 
quantitative analysis. 

To summarize, our system of equations reduces to 
Eqs. ©-Q and The "control parameter" in sub- 

sequent analysis will turn out to be 



w 



A 

r 



(12) 



which, when fixed, leaves the depression factor r as an 
overall control of rate of plasticity. 
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III. MEAN-FIELD APPROXIMATION 
A. Single neuron response function 

To formulate a mean-field approximation, we assume 
that the firing rates of the neurons are given by the same 
mean-field value A and that all the synapses have the 
same synaptic weight w. The total synaptic conductance 
G of a singled-out postsynaptic neuron is a function of 
time and jumps whenever there is a firing of a presynaptic 
neuron and subsequently decays with the time constant 
Tm as in Eq. ([3|) . Consider the limit of large number K of 
afferent synapses per neuron while keeping the product 
luK constant. Then, the total frequency of all presynap- 
tic firings is KX —^ oo. In this limit, the jump sizes of 
G, being proportional to w, approach while the rate of 
making the jumps diverges. For fixed A and wK, we may 
thus ignore fluctuations and consider G as a constant over 
time. 

For a constant total synaptic conductance G, the fir- 
ing frequency of an integrate-and-fire neuron can be 
solved exactly by setting the initial membrane potential 
to (0) — Vi, the reset value, and finding the time r 
it takes to reach the firing threshold V {t) = Vth- The 
solution of Eq. ^ is given by 



G + 1 



■In 



V (G) - -j/th 

V (G) - Vr 



when V (G) > Vth, where 
1/(G)e 



Vq + GTZ 
G + 1 



(13) 



(14) 



is the resting value for the membrane potential under the 
constant conductance G. For V < Vth, we have t — oo 
since the neuron never fires. In addition to firings due to 
crossing the threshold, the period t can be cut short by 
Poisson noise of frequency Aat as the membrane poten- 
tial sweeps from Vi- to Vth- Given Poisson statistics, the 
average interval r between all firings (either threshold- 
crossing or noise-driven) is given by 



dtXNte 



X-^Ul-e-^«^), 



(15) 



and the firing frequency A of the neurons is given by 



A(G) = Aat 



V (G) - Vth 
V{G)-V, , 



(16) 



for V > Vth, and by A = Aat otherwise. (See Fig.[T]) The 
function A (G) obtained in (|16p characterizes the mean- 
field response of the integrate-and-fire network model in 
the current study j45j . but similar functions can be ob- 
tained analytically or numerically for different neuron 
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Figure 1: Firing frequency A (G) of an integrate-and-fire neu- 
ron under constant total synaptic conductance G given by 
Eq. (I16|) . Also shown are mean-field active synaptic trans- 
mitter fractions 1^ (A) (for K — 31, w — 0.01) driven by a 
Poisson spike train of frequency A given by (|17p . The dashed 
line represents a numerically calculated correction to Y (A) 
due to the periodic firing of the neuron when driven by con- 
stant total synaptic conductance G. The parameters are listed 
in Table |T] The arrow indicates that the curve shifts to the 
right for increasing w. The dot-dashed line is the neuron re- 
sponse function numerically calculated for the Morris-Lecar 
model used in [sj with a slightly higher background current. 



Table I: Values of parameters used in calculations 



integrate & fire 



TUM model 



resting potential Vb 

leak time 
firing threshold 

reset potential K 
reversal potential TZ: 



-55 mV 
20 ms 
-54 mV 
-80 mV 
mV 



decay time td- 20 ms 

recovery time tr: 200 ms 

release fraction u: 0.5 

noise frequency Ajv: 1 Hz 

plasticity rate r: 0.01 



models, such as, the Hodgkin-Huxley, Morris-Lecar, or 
FitzHugh-Nagumo models (see, e.g., [HI [11]). (An ex- 
ample of a numerically calculated response function for 
a Morris-Lecar neuron is included in Fig. [T]) 



B. Synapse response function 

To obtain the mean-field active transmitter fraction Y 
of a neuron, we further approximate the firing of the neu- 
rons as having Poisson statistics described by the mean- 
field firing rate A [l3|. Keeping Y and Z as the inde- 
pendent variables for the TUM model, the stochastic dy- 
namics ^ can be averaged over an ensemble of Poisson 
spike trains with (S) = X for the presynaptic neuron to 
yield the functions 



Y{X) = 



UTjjX 



1 + u{td + tr) X 



(17) 



and Z (A) = —Y (A) for the average active {Y) and inac- 
tive (Z) transmitter fractions. The function Y (A) char- 
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synaptic weight w 

Figure 2: Mean-field firing frequency A as a function of mean- 
field synaptic weight w determined by the intersections of neu- 
ron response function and synapse response function as plot- 
ted in Fig. [T]for K — 31. The two arrows show the region of 
hysteresis where two stable fixed points, the noise-dominating 
lower fixed point and the persistently active upper fixed point, 
coexist. The symbols will be used in connection with Figs. |4] 
and El 



acterizes the mean- field response of a synapse [45]. The 
value of the active transmitter fraction for a mean-field 
network can be obtained through the relation 



Y = Y{X). (18) 



C. Self-consistency condition 

The total synaptic conductance of a neuron in a mean- 
field network with K afferent synapses per neuron is given 
by 



G = KwY = KwY (A) , 



(19) 



which can be substituted into Eq. (|16p to complete the 
self-consistency condition closing the set of Eqs. 
p7)) . and The curve for Eq. ^ is also shown 

in Fig. [1] and the intersections with the neuron response 
function A (G) represent fixed points of the network ac- 
tivities, A, given the mean-field synaptic weight w. The 
function A (w) is shown in Fig. [2l The mean-field the- 
ory predicts a first-order phase transition and hysteresis 
as the mean-field synaptic weight w is varied. At low 
synaptic weight, the activities of the network are dom- 
inated by the external noise. As the synaptic weight is 
increased, a pair of new fixed points emerge at higher 
firing frequencies. Among the three fixed points, only 
the upper and the (noise-dominated) lower fixed points 
are stable. As the synaptic weight is increased further, 
the lower fixed point eventually is annihilated by the un- 
stable (middle) fixed point leaving only the upper fixed 
point representing higher firing activities. 




Figure 3: We consider the dynamics of a single synapse of 
synaptic weight w in a mean-field environment. 



IV. RANDOM WALK IN SYNAPTIC WEIGHT 
SPACE 

Given the mean-field firing frequency A for a fixed 
mean synaptic weight w, motivated by the Bethe-Peierls 
approximation [25|, we consider a single synapse within 
such a mean-field environment as illustrated in Fig. [31 
The time-dependence of the synaptic weight is governed 
by Eq. ^ and the simplification embodied in Eq. (|lip . 
With the assumption of Poisson spike trains produced 
by both the presynaptic and postsynaptic neurons, the 
dynamics of the synaptic weight maps onto an asym- 
metric random walk (jump process) in w-space with w- 
dependent jump rates and jump sizes [37| . From Eqs. ([H) 
and (|lip , the synaptic weight increases by A A2 = A I2 
whenever the postsynaptic neuron fires and decreases by 
rwBi = rwYi whenever the presynaptic neuron fires. As 
the presynaptic neuron is driven by the mean-field envi- 
ronment, its firing frequency A2 is given by the mean-field 
firing frequency A. On the other hand, the postsynap- 
tic neuron is partly driven by the synapse in question 
(Fig. |3l), thus its firing frequency Ai is w dependent. As 
calculated from the neuron response function A (G), this 
firing rate is given by 



Ai = X{[{K -l)w + w] Y) 



(20) 



where, even though the synaptic weight w is singled out 
from among the K afferent synapses, all of the presy- 
naptic neurons are from the mean-field environment (see 
Fig. [3]) with the active transmitter fraction Y. As for the 
step sizes, since A2 depends on the firing of the presy- 
naptic neuron, its value is given by the mean-field value 
Y. On the other hand, Bi is determined by the firing of 
postsynaptic neuron; therefore its value is also w depen- 
dent and given by Bi = y(Ai). As mentioned earlier, 
we set A ^ B ~ Y to reduce the number of variables and 
computational intensiveness. In general, the system can 
have different characteristic functions A (A) and B (A) for 
potentiation and depression. 



A. Mean- field synaptic weight 

While the "random walk" of synaptic weight is un- 
der the influence of the mean-fleld environment, a self- 
consistency condition requires the average value of the 
singled-out synaptic weight from the random walk pro- 
cess to be the same as the assumed mean-field value 
w. Assuming a normalized stationary synaptic-weight 
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Figure 4: Average synaptic weight from the simulation of 
a single synapse between two Poisson neurons as shown 
schematically in Fig. |3] versus the mean-field synaptic weight 
w for w* = 0.02. The mean-field firing frequencies A used in 
the single synapse simulations are given by Fig. [2] and the 
results for both the upper and lower fixed points (marked by 
corresponding symbols in Fig. [2| are plotted. 



distribution from the random walk, Ps (w), which is w- 
depcndcnt, the self-consistency requires 



— (w) — / wPs {w; w) dw. 



(21) 



The stationary synaptic weight distribution, following 
the stochastic dynamics (fTTI) . can not be easily deter- 
mined even with the assumption of Poisson spike trains. 
However, straightforward numerical simulations of a sys- 
tem of a single synapse and two Poisson neurons can be 
used to obtain the stationary distribution Pg (w; w) to 
any reasonable precision. Figure [¥] shows typical plots 
of the average synaptic weight from the random walk 
process as a function of the mean-field synaptic weight 
w. Notice that the symbols are lower (higher) than the 
dashed line on the left (right) meaning that a deviation of 
w from w* results in a smaller deviation of (w) from w* . 
Hence, the w = w* is a. stable solution for the condition 



B. Fokker— Planck equation 

Since the synaptic weight dynamics ^ [or ([TT]) ] only 
depend on the current weight of the synapse (a Markov 
process) , we can approximate the dynamics of the weight 
distribution for small jumps with a master equation ig- 
noring higher-order moments in the Kramers-Moyal ex- 
pansion (see, e.g., [37]). A Fokkcr-Planck equation for- 
malism has been applied to populations of neurons in 
previous studies of the distributions in membrane poten- 
tial [13, Ell, H^l ■ Here we apply this type of formalism to 
a population of plastic synapses similar to the work by 
Rubin, Lee, and Sompolinsky p7j . 



Equation ([8]) or its simplification ((TT|) can be viewed 
as a stochastic equation governing the jumps in synaptic 
weight space upon arrival of spikes and which allows us 
to calculate the drift 



lim 

At^O 



{w (At) - w) 



w(0)—w 



At 

rwY{\i)\, 



(22) 



where Ai is given by Eq. ([20]) . One also has for the dif- 
fusion coefhcient (arising from the second moment of the 
synaptic weight changes) 



D 



— lim 

2 At^o 



w (At) 



w(0)—w 



At 



= ^A^Y^X, + ^r'w^Y{X,f~X. (23) 

The dynamics of the synaptic weight distribution P (w) is 
then approximated for small jumps by the Fokker-Planck 
equation 

_|_[„P(^)] + |!_[Z?P(^)] (24) 

when higher-order moments in the Kramers-Moyal ex- 
pansion are ignored. Consider the stationary state P^ (w) 
of the weight distribution. The Fokker-Planck equation 
can be integrated once to yield 



dP{w) 
dt 



VPs {w) = ^ [DPs [w) 
whicli is formally solved by 



const.. 



(25) 



(26) 



The peak positions of tlic distribution (|26p . w, (if they 
exist) are given by v{'w) — D' (w). Since the diffusion 
coefficient (|23l) is one-order higher in the plasticity rate 
r compared to the drift (|22p. under the small jump ap- 
proximation, the peak positions are given approximately 
by the zeros of v. When w = w* = A/r, we can ver- 
ify that the drift is zero and thus the distribution 
Ps (w) should peak at w = w. We can expand v and D 
in powers oi w — w and keep only the first order terms. 
The exponent of Eq. ((26l) becomes 



dw'- 



rYw* 



{w* — wY 
2rYw*'^ 



const.. 



(27) 



which results in a Gaussian distribution for Ps (w). The 
distribution is very narrow; the width of the distribution 
is given by 

Aw ~ \/r?w* (28) 

which corresponds to only about 1% oi w* for the pa- 
rameters we have considered in Table |T] above. Simula- 
tions for a fully connected network using the same model 
and parameters yield a width about five times larger but 
with the same scaling power in r in the noise-dominated 
regime 53] . We will return to this point below. 



8 



C. Stability analysis within mean- field theory 

The condition v (w) = for the peak position w of the 
synaptic-weight distribution gives rise to 



AyAi = rwY 



(Ai) a, 



(29) 



where Ai = Ai (w). Using the mean-field response func- 
tion pT]) for the synapses, we get 



w = w 



1 -I- M {td + tr) Ai 
1 + u{td + tr) a 



(30) 



It is straightforward to verify that when the mean-field 
synaptic strength is given by w ^ w* , the peak position 
is given hy w = w, providing a self-consistent solution. 
To determine the stability of the solution, we start with 
a small deviation Sw of w from w* and find the resultant 
deviation Sw of the peak position w from w. Assum- 
ing the neuron response function A (G) is smooth around 
the mean-field total conductance G = KwY, the firing 
frequency of the postsynaptic neuron can be expanded 



Ai w A + A' (G) dwY 



(31) 



where A' is the first derivative of the mean-field response 
function A(G). The mean- field response function Y {X) 
for the active transmitter can also be expanded as 



Y{Xi) = Y[l + Y' (A) A' (G) Sw] 
The condition v (w) ~ results in 



(32) 



w-w* = bw^\ w*\' (G) 



-l}Sw. (33) 



For the w — w* solution to be stable, 5w and Sw should 
be of opposite signs so that the perturbation can be 
damped. Thus, we need 



w*\' (G) 



Y 



y'(A) 



< 1, 



or, for given total synaptic conductance G, 



v*<wt (G)^ A' (G) 



Y 
J 



r'(A) 



(34) 



(35) 



For the integrate-and-fire model we have considered, in 
the noise-dominated regime A = Xn and A' — 0, so the 
condition (|35|) is satisfied. In the large G limit, the firing 
frequency A increases linearly with G while the trans- 
mitter fraction Y for the TUM model saturates leading 
to a linear w^ (G) . For self-consistency, the membrane 
conductance is given by G = Kw*Y for the w = w* 
solution, and it is thus left to the constant K (represent- 
ing the number of connections per neuron) to determine 
whether the solution will remain stable. For the range of 



the parameters we have considered, the solution w = w* 
remains stable for large w* as long as > 1. However, 
it is straightforward to verify from Eq. ([35]) that iJu — w* 
can not remain a stable solution at large w* for mod- 
els in which, for example, A (G) saturates for large G or 
Y (A) increases linearly (when it no longer represents a 
fraction) with A for large A. Under these situations, this 
analysis suggests that one may find runaways or pile-ups 
of the synaptic weights in the system. 



V. CORRECTION TO THE MEAN-FIELD 
THEORY 

In the mean-field considerations, we characterized all 
neurons with a single firing frequency A, all synapses with 
a constant active transmitter fraction Y, and the "net- 
work" by a constant strength w. This is an over sim- 
plification even within the mean-field approach. We re- 
tain the characterization of the environment by a single 
w. However, a neuron within a network experiences, in- 
stead of a constant total synaptic conductance, the bom- 
bardment of synaptic conductance pulses issuing from 
the spike trains of corresponding presynaptic neurons. 
This shot-noise-like fiuctuation is most prominent when 
the amplitude of the pulses is comparable with the mean 
of the total synaptic conductance. When the number K 
of afferent synapses per neuron and the total frequency 
Atotai of presynaptic events are small, the mean level of 
the total synaptic conductance is comparable with the 
amplitude of each individual pulse. Consequently, the 
assumption of a constant total synaptic conductance will 
be inadequate. 

In the mean-field approximation above, we also as- 
sumed that the neurons fire following Poisson statistics. 
However, when integrate-and-fire neurons are driven by 
a constant total synaptic conductance (as in the limit of 
large K and Atotai), the time between threshold crossings 
will also be constant, and the resulting spike trains of the 
neuron will be periodic. Combined with the Poisson ex- 
ternal noise, the intervals between firings will follow a 
Poisson distribution only up to the period of threshold 
crossings, where there will be a (5-function representing 
the periodic firings. Nonetheless, Poisson and periodic 
firings result in the same mean level of active transmit- 
ter fraction in the limits of low and high firing frequen- 
cies, while the corrections in the intermediate frequency 
regime remains insignificant for the TUM model that we 
considered (see the dashed line in Fig. [T|). We will thus 
retain the assumption of Poisson firing for the following 
discussion. 



Effect of fluctuation in total synaptic 
conductance 



Apart from making an analytically tractable system, a 
mean-field approach does not necessarily require approx- 
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Figure 5: Extended space of Fig. [T] when the fluctuation 
in synaptic conductance G is considered (represented by the 
jump size J). The extended neuron response function (|38p . 
the synapse response function Y (A) = G/ (Kw), and the re- 
lation (|39l) determine the stationary states of a mean-fleld 
network (marked by a circle on the surfaces). 




' ' ' ' I I I I I ' ' ' ' I I I I I ' ' ' ' I 

synaptic weight w 



Figure 6: Numerical mean-field phase diagram (firing fre- 
quency A versus synaptic weight w) for various numbers K (as 
labeled in legend) of afferent synapses per neuron under fluc- 
tuating total synaptic conductance (drawn lines) compared 
with the analytical results when the total synaptic conduc- 
tance G = G is assumed to be constant (dashed line). 



imating the total synaptic conductance as a constant in 
time. Here, we consider the total synaptic conductance 
G as a time dependent function and model the dynamics 
of G approximately with 



~r~ — h J iS'total, 

at td 



(36) 



where S'totai is a Poisson spike train with the frequency 
Atotai accounting for any presynaptic events, and the 
jump size J is a constant modeling the increase of the 
total synaptic conductance due to each spike in S'totai- 
(In the spirit of mean-field theory, we simplify the dy- 
namics of G by assuming a constant jump size J. In gen- 
eral, the jump size of the total synaptic conductance for 
each presynaptic event is a random variable depending 
on the active transmitter fraction of the firing presynap- 
tic neuron.) The average total synaptic conductance G 
is related to frequency Atotai and jump size J by 



G = AtotalTuJ, 



(37) 



and the ensemble of G (t) can be characterized by two 
independent parameters from the set Atotai, J, and G. 
Here, we choose G and J to characterize the ensem- 
ble, and the resulting ensemble-averaged firing frequency 
Aout {G, J) of the post-synaptic neuron is a generalization 
of the response function A (G) considered in Section Hill 
To evaluate Aout (G, j) for given values of G and J, 
we simulate the dynamics p6p with a Poisson spike train 
of frequency Atotai inferred from Eq. (j37p to obtain the 
time-dependent total synaptic conductance G (t) , which 
we use in the dynamics of the membrane potential ([3|: 
the measured average firing frequency of the postsynaptic 
neuron gives us the value of Aout (G, j) . The process is 
repeated to generate the surface 



A = Ao 



(G,J) 



(38) 



as plotted in Fig. \5\ compare to the line representing 
A (G) as plotted in Fig. [T] The jump size J is an estima- 
tor of the fluctuations in the total synaptic conductance 
G (t), and in the J = limit, the two-parameter response 
function Aout (G, A) reduces to the single parameter re- 
sponse function Aout (G, O) = A (G) . 

As we are considering a mean-field network with K 
afferent synapses per neuron, the total frequency of 
presynaptic events for a postsynaptic neuron is given by 
Atotai = KX, where A is the mean-field firing frequency 
for any neuron. Combining with the relation (|37|) . this 
gives us the condition 



G 



KtdJ 



A, 



(39) 



as plotted in Fig. [5l forming a planar surface in the 
logarithmic-scale G~J-X space. 

Whatever the value of J, the synaptic response func- 
tion Y (a) determines the mean total synaptic conduc- 
tance G = KwY (a) when A is given for the mean-field 
network. The inverse 



\Kw J 



(40) 



also plotted in Fig. [S] with Eq. , represents the con- 
vex surface invariant along the J-axis. Along with the 
surfaces and ([M)) . the intersections of the three sur- 
faces determine the fixed points of the system (marked 
with a circle in Fig. [5]). While the neuron response func- 
tion Aout (G, j) is evaluated numerically, we determine 
the fixed points in the mean-field firing frequency nu- 
merically from the intersections of the three surfaces in 
Fig. [5l varying the mean synaptic weight w for a given 
number of afferent synapses K per neuron. The resulting 
phase diagram is plotted in Fig. [S] for various K values. 
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The hysteresis vanishes at K = Kc ~ 10, and the transi- 
tion from noise dominance to persistent activity becomes 
continuous. In our case, the fluctuation comes from the 
fact that the number of afferent synapses is finite. We 
expect that, in a real biological network, other sources 
of fluctuations can play a similar role to smooth out the 
first-order phase transition separating a noise dominated 
regime from persistent activity as predicted by the ana- 
lytical mean- field theory of Section UTTl 

B. Correction to synaptic weight distribution 

Besides smoothing out the first-order phase transition 
in a mean-field network, the fluctuation in total synaptic 
conductance can also influence the stationary distribu- 
tion of synaptic weights. We have noted that the width 
of the synaptic weight distribution predicted by the sim- 
plest mean-field approach (about 1% of w*) is too small 
in comparison with results of simulations on fully con- 
nected plastic integrate-and-fire networks (about 5% of 
w* [43II). One of the approximations we made in arriving 
at Eq. (|27|) is to replace the time-dependent Y (or the 
fractions A and B that it represents) with its average 
value. Such an approximation can be improved by con- 
sidering the second moment of Y in deriving the diffusion 
coefficient D in Eq. (|23p . We expect the fluctuation in 
Y to be important when the mean-field firing frequency 
is low. In this limit the time-dependent Y (t) is a sum of 
pulses 

y(t)=7.^0(i-t,)e-(*-*')/^-, (41) 

i 

where 9 (t) is a step function with 9 {t > 0) — 1 and 
6 {t < 0) — 0, {ti} are the firing times of the neuron, and 
mean (Y) — utdX. When the overlap of the pulses can be 
ignored, for example, in the noise dominated regime, the 
mean square of Y is given by {Y'^) = u^T£)A/2, which 
should replace the squares of Y in the diffusion coeffi- 
cient ([23t . This leads to a width of the synaptic weight 
distribution of 



Aw ~ y r^^Y^w* = ^/ur/2w* (42) 

instead of Eq. . For the values of the parameters we 
have considered, this represents about 5% of w* and is 
similar to the simulation results from a fully connected 
network ^43| in the noise-dominated regime. We have 
verified the analytical result (j42p through numerical sim- 
ulations of a random walk process for the synaptic weight 
described by the dynamics pT|) . As shown in Fig. [71 the 
scaling in plasticity rate r continues to be described by 
Aw ^ r" with a ~ 1/2 in all cases, while the amplitudes 
match well for the noise-dominated fixed point but show 
small deviations when the system is persistently active. 
In general, we can expect fluctuations to play an impor- 
tant role in determining the synaptic weight distribution 




ID 



10"^ I ■ ^ 

10"-^ 10"^ 10"' 

plasticity rate r 

Figure 7: Numerical results (symbols) for the width of 
synaptic-weight distribution, /S.w, as a function of plastic- 
ity rate, r, in a limited system of a single synapse between 
two neurons with Poisson spike trains following the dynamics 
pip and operating at three of the fixed points predicted by 
the mean-field theory (as marked by corresponding symbols 
in Fig. [21). The results are compared with the analytical pre- 
dictions (lines of corresponding colors) from Eq. (|42p (drawn 
lines) and Eq. (1^ (dashed line). 



of a plastic network. Here we have shown that a signifi- 
cant correction can be accounted for by considering the 
fluctuation in the neural transmission for a limited sys- 
tem of a single synapse between two neurons with Poisson 
spike trains. Even without an analytical solution, such a 
system can be easily simulated to any desirable accuracy 
to obtain the stationary synaptic weight distribution as 
well as the average synaptic weight (w), which provides 
the self-consistency condition pi]) for the mean-field the- 
ory. 

Here we have preserved the basic structure of the 
mean-field approach under the additional consideration 
of fluctuations in synaptic conductance by extending the 
neuron and synapse response functions to functions of 
two variables. The additional degree of freedom, the av- 
erage jump size J of active transmitter fractions, was 
fixed by the condition ([551) coming from the nature of 
such a fluctuation and does not enter the final result (|42p 
explicitly. 



VI. SUMMARY AND CONCLUSION 

Models of plastic neural networks can generally be bro- 
ken down into three parts: (i) the modeling of the dynam- 
ics of the neuron state (represented by the membrane po- 
tential), (ii) the modeling of synaptic transmission, and 
(iii) the plasticity rules. We have applied a mean-field 
theory that follows this basic structure. In the mean- 
field framework of the current study, we reduce the dy- 
namics of the neuron state to a response function A (G) 
representing the mean firing frequency of a neuron when 
it is driven by a constant total synaptic conductance G. 
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Similarly, the dynamics of synaptic transmission is also 
reduced to a characteristic response function Y (A) rep- 
resenting the mean fraction of active neural transmitter 
given the Poisson firing frequency A of the presynaptic 
neuron. With the mean-field synaptic weight w of a 
network with K afferent synapses per neuron and the 
mean- field firing frequency A, the average total synaptic 
conductance is given by G = KwY (A). This allows one 
to plot the synapse response function Y (A) along with 
the neuron response function A (G) (see Fig. [T|) to deter- 
mine self-consistently the mean-field firing frequency A 
as a function of w and K. 

The mean-field firing frequency A and the mean synap- 
tic weight H) characterize the mean-field network com- 
pletely. We then use this mean-field network as an en- 
vironment and investigate the dynamics of spike-timing- 
dependent plasticity ([8]) of a single synapse in such an 
environment. Assuming Poisson statistics for the spike 
trains, the dynamics ^ can be viewed as describing a 
type of random walk in synaptic weight space, where the 
frequency of potentiation and the jump size of depression 
are dependent on the weight of this synapse. The distri- 
bution of synaptic weights can be calculated numerically 
for the stationary state using straightforward simulations 
of the single synapse system. Under the approximation 
of small jumps, the stationary distribution can be calcu- 
lated analytically from the Fokker-Planck equation aris- 
ing from the master equation. The self-consistency of the 
mean-field approach is completed by requiring that the 
mean of the synaptic weight distribution (w) reproduces 
the mean- field synaptic weight w of the environment. 

We have considered a network of integrate-and-fire 
neurons 13511 coupled through synapses with TUM dy- 
namics [27| within the mean-field framework outlined 
above. We chose integrate-and-fire neurons for their sim- 
plicity and wide use. We chose the TUM model of neu- 
ral transmission for its features and flexibility as noted 
above. On the simplest level for a static network, our 
mean-field approach amounts to finding intersects of the 
neuron response function and synapse response function, 
each of which can be calculated separately from the par- 
ticular neuron model and synapse model selected. Our 
choice of model happens to allow us to find analytical 
forms of the mean-field response functions for the neu- 
rons and synapses, and predict a first-order phase tran- 
sition from a noise-dominated regime to a regime of per- 
sistent activity as the mean-field synaptic weight is in- 
creased. (However, see below.) In general, these response 
functions can be computed numerically in a straightfor- 
ward fashion for a variety of neuron and synapse models 
regardless of the number of empirical parameters they 
might carry. In Fig. [H we showed the results of such cal- 
culation for the Morris-Lecar neuron used in [s^ which 
was defined by more than twenty parameters. 

For the dynamics of the synaptic weight, for speci- 
ficity, we follow the spike-timing-dependent plast icity 
rules proposed by van Rossum, Bi, and Turrigiano [22l |. 
with additive potentiation and multiplicative depression. 



The Fokker-Planck analysis of the corresponding ran- 
dom walk process predicts a narrow Gaussian distribu- 
tion for the synaptic weight centered around w*, the con- 
trol parameter entering the plasticity rules. [See, 
e.g., Eq. (|lip .] We apply a small perturbation to the 
{w) — w — w* solution, analyze its stability analytically 
and numerically, and find it to be stable for any number 
of afferent synapse per neuron K > I for the model and 
parameter ranges we considered. However, the analyti- 
cal expression psp for the stability, which is generic for 
any neuron and synapse response functions, also suggests 
that for models in which the neuron firing rate A (G) satu- 
rates for large conductance G (for example, models with 
a refractory period for the firing of neurons) or where 
the synapse response function Y (A) does not saturate 
for large A, (for example, when the effect of presynaptic 
firings is additive and Y no longer represents a fraction) , 
the fixed point id = w* can not remain stable for large 
w* . It is then possible to have runaway or pile- up in the 
resulting synaptic weight distribution. 

For a network of finite K and low overall frequency 
Atotai of presynaptic events, the "shot noise" due to the 
discrete nature of presynaptic spikes is not negligible, 
and one at least needs to expand the description of the 
mean-field response of a neuron to include temporal fluc- 
tuations. We have done this approximately via a two- 
variable function, e.g., Aout (G,J), where J, represent- 
ing fluctuations, describes the size of the jumps in a 
neuron's total synaptic conductance for each presynap- 
tic event. We evaluate this two-variable response func- 
tion numerically for a single neuron, modeling its total 
synaptic conductance with a simple stochastic jump-and- 
decay process l\'S6\i . The results suggest the disappear- 
ance of the first-order transition when the number, K, of 
afferent synapses per neuron is less than a critical value 
Kc ~ 10. In realistic situations fluctuations can be ex- 
pected to smear out the first-order transition. When a 
corresponding fluctuation is considered for the variable 
Y governing jump sizes in the plasticity rules pip , the 
Fokker-Planck approach predicts a broader Gaussian dis- 
tribution which has a width similar to the observed width 
in full network simulations [i^ . The results of extensive 
simulations on a fully connected plastic network will be 
published elsewhere. 

The inclusion of fluctuations at some level, such as 
within our extended mean-field theory (see Sec. |V|, hints 
at a possible critical state of the system at the endpoint 
of a first-order transition line in an alog y with the vapor 
pressure curve of a fluid, see, e.g., [2^. While critical- 
ity in neural avalanches has been observed by Beggs and 
Plenz 31], within the extended mean- field analysis em- 
ployed in the current study the system does not appear 
to organize into such a critical state without a requisite 
tuning of the fluctuation-amplitude and perhaps plastic- 
ity parameters, e.g., w* . This is in contrast with the 
suggestion that such criticality should be self-organized 
[Sll . HJ. It thus will be of great interest to find miss- 
ing elements, possibly in the plasticity rules, that could 
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dynamically push the network close to a critical state. 
Along this direction, one can expect a variety of model 
candidates, which can be easily subjected to the type 
of mean-field scheme outlined in the current study to 
provide additional qualitative and semi-quantitative in- 
sight into their plausibility. We also note in passing that 
"all-to-all" network simulations using the plasticity rules 
and neuron modeling described in this paper have not 
revealed evidence for a self-organized critical state. In 
that regard, sparse networks would appear to be better 
suited, but mean-field analysis such as presented here 
will mainly have only qualitative use. Our simulations of 
integrate-and-fire and other neuronal networks with the 
plasticity rules used here will be presented elsewhere. 

As noted in the Introduction, the current study cannot 
begin to address how brains function or form. As a real 
brain is never uniform or stationary, we do not expect 
the model systems presented here to address or repro- 
duce its dynamics. However, there are cultured networks 



consisting of hundreds of neurons with virtually "all-to- 
all" interactions [H, [H, [s^, [U, [12] . The studies of these 
networks in vitro is an important steppingstone towards 
the understanding of more complicated networks. The 
dynamics of these cultured network are on a scale very 
much accessible to current neural network modeling, such 
as the one presented here, and simulation, as noted, to 
be presented elsewhere. 
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